STA 364 Time Series Applications

Author

Affiliation

Ch8. Exponential smoothing

?var:college.name
?var:course.number ?var:course.semester

Exponential smoothing

Historical perspective

• Developed in the 1950s and 1960s as methods (algorithms) to produce point forecasts.
• Combine a “level”, “trend” (slope) and “seasonal” component to describe a time series.
• The rate of change of the components are controlled by “smoothing parameters”: \(\alpha\), \(\beta\) and \(\gamma\) respectively.
• Need to choose best values for the smoothing parameters (and initial states).
• Equivalent ETS state space models developed in the 1990s and 2000s.

Big idea: control the rate of change

\(\alpha\) controls the flexibility of the level (Level - the average value for a specific time period)

• If \(\alpha = 0\), the level never updates (mean)
• If \(\alpha = 1\), the level updates completely (naive)

\(\beta\) controls the flexibility of the trend

• If \(\beta = 0\), the trend is linear
• If \(\beta = 1\), the trend changes suddenly every observation

\(\gamma\) controls the flexibility of the seasonality

• If \(\gamma = 0\), the seasonality is fixed (seasonal means)
• If \(\gamma = 1\), the seasonality updates completely (seasonal naive)

A model for levels, trends, and seasonalities

We want a model that captures the level (\(\ell_t\)), trend (\(b_t\)) and seasonality (\(s_t\)).

ETS models

Additive ("A") or multiplicative ("M")

None ("N"), additive ("A"), multiplicative ("M"), or damped ("Ad" or "Md").

None ("N"), additive ("A") or multiplicative ("M")

Simple exponential smoothing

Simple methods

Time series \(y_1,y_2,\dots,y_T\).

• Want something in between these methods.
• Most recent data should have more weight.

Simple Exponential Smoothing

Simple Exponential Smoothing

Simple Exponential Smoothing

• \(\ell_t\) is the level (or the smoothed value) of the series at time t.
• \(\pred{y}{t+1}{t} = \alpha y_t + (1-\alpha) \pred{y}{t}{t-1}\)

Iterate to get exponentially weighted moving average form.

Optimising smoothing parameters

• Need to choose best values for \(\alpha\) and \(\ell_0\).
• Similarly to regression, choose optimal parameters by minimising SSE: \[ \text{SSE}=\sum_{t=1}^T(y_t - \pred{y}{t}{t-1})^2. \]
• Unlike regression there is no closed form solution — use numerical optimization.

• For Algerian Exports example:
  • \(\hat\alpha = 0.8400\)
  • \(\hat\ell_0 = 39.54\)

Simple Exponential Smoothing

Models and methods

Methods

• Algorithms that return point forecasts.

Models

• Generate same point forecasts but can also generate forecast distributions.
• A stochastic (or random) data generating process that can generate an entire forecast distribution.
• Allow for “proper” model selection.

ETS(A,N,N): SES with additive errors

Forecast error: \(e_t = y_t - \pred{y}{t}{t-1} = y_t - \ell_{t-1}\).

Specify probability distribution for \(e_t\), we assume \(e_t = \varepsilon_t\sim\text{NID}(0,\sigma^2)\).

ETS(A,N,N): SES with additive errors

where \(\varepsilon_t\sim\text{NID}(0,\sigma^2)\).

• “innovations” or “single source of error” because equations have the same error process, \(\varepsilon_t\).
• Measurement equation: relationship between observations and states.
• State equation(s): evolution of the state(s) through time.

ETS(M,N,N): SES with multiplicative errors.

• Specify relative errors \(\varepsilon_t=\frac{y_t-\pred{y}{t}{t-1}}{\pred{y}{t}{t-1}}\sim \text{NID}(0,\sigma^2)\)
• Substituting \(\pred{y}{t}{t-1}=\ell_{t-1}\) gives:
  • \(y_t = \ell_{t-1}+\ell_{t-1}\varepsilon_t\)
  • \(e_t = y_t - \pred{y}{t}{t-1} = \ell_{t-1}\varepsilon_t\)

• Models with additive and multiplicative errors with the same parameters generate the same point forecasts but different prediction intervals.

ETS(A,N,N): Specifying the model

ETS(y ~ error("A") + trend("N") + season("N"))

By default, an optimal value for \(\alpha\) and \(\ell_0\) is used.

\(\alpha\) can be chosen manually in trend().

trend("N", alpha = 0.5)
trend("N", alpha_range = c(0.2, 0.8))

Example: Algerian Exports

algeria_economy <- global_economy %>%
  filter(Country == "Algeria")
fit <- algeria_economy %>%
  model(ANN = ETS(Exports ~ error("A") + trend("N") + season("N")))
report(fit)

Series: Exports 
Model: ETS(A,N,N) 
  Smoothing parameters:
    alpha = 0.8399875 

  Initial states:
   l[0]
 39.539

  sigma^2:  35.6301

     AIC     AICc      BIC 
446.7154 447.1599 452.8968 

Example: Algerian Exports

components(fit) %>% autoplot()

Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_line()`).

Example: Algerian Exports

components(fit) %>%
  left_join(fitted(fit), by = c("Country", ".model", "Year"))

# A dable: 59 x 7 [1Y]
# Key:     Country, .model [1]
# :        Exports = lag(level, 1) + remainder
   Country .model  Year Exports level remainder .fitted
   <fct>   <chr>  <dbl>   <dbl> <dbl>     <dbl>   <dbl>
 1 Algeria ANN     1959    NA    39.5    NA        NA  
 2 Algeria ANN     1960    39.0  39.1    -0.496    39.5
 3 Algeria ANN     1961    46.2  45.1     7.12     39.1
 4 Algeria ANN     1962    19.8  23.8   -25.3      45.1
 5 Algeria ANN     1963    24.7  24.6     0.841    23.8
 6 Algeria ANN     1964    25.1  25.0     0.534    24.6
 7 Algeria ANN     1965    22.6  23.0    -2.39     25.0
 8 Algeria ANN     1966    26.0  25.5     3.00     23.0
 9 Algeria ANN     1967    23.4  23.8    -2.07     25.5
10 Algeria ANN     1968    23.1  23.2    -0.630    23.8
# ℹ 49 more rows

Example: Algerian Exports

fit %>%
  forecast(h = 5) %>%
  autoplot(algeria_economy) +
  labs(y = "% of GDP", title = "Exports: Algeria")

Models with trend

Holt’s linear trend

• Two smoothing parameters \(\alpha\) and \(\beta^*\) (\(0\le\alpha,\beta^*\le1\)).
• \(\ell_t\) level: weighted average between \(y_t\) and one-step ahead forecast for time \(t\), \((\ell_{t-1} + b_{t-1}=\pred{y}{t}{t-1})\)
• \(b_t\) slope: weighted average of \((\ell_{t} - \ell_{t-1})\) and \(b_{t-1}\), current and previous estimate of slope.
• Choose \(\alpha, \beta^*, \ell_0, b_0\) to minimise SSE.

ETS(A,A,N)

Holt’s linear method with additive errors.

• Assume \(\varepsilon_t=y_t-\ell_{t-1}-b_{t-1} \sim \text{NID}(0,\sigma^2)\).
• Substituting into the error correction equations for Holt’s linear method \[\begin{align*} y_t&=\ell_{t-1}+b_{t-1}+\varepsilon_t\\ \ell_t&=\ell_{t-1}+b_{t-1}+\alpha \varepsilon_t\\ b_t&=b_{t-1}+\alpha\beta^* \varepsilon_t \end{align*}\]
• For simplicity, set \(\beta=\alpha \beta^*\).

Exponential smoothing: trend/slope

ETS(M,A,N)

Holt’s linear method with multiplicative errors.

• Assume \(\varepsilon_t=\frac{y_t-(\ell_{t-1}+b_{t-1})}{(\ell_{t-1}+b_{t-1})}\)
• Following a similar approach as above, the innovations state space model underlying Holt’s linear method with multiplicative errors is specified as \[\begin{align*} y_t&=(\ell_{t-1}+b_{t-1})(1+\varepsilon_t)\\ \ell_t&=(\ell_{t-1}+b_{t-1})(1+\alpha \varepsilon_t)\\ b_t&=b_{t-1}+\beta(\ell_{t-1}+b_{t-1}) \varepsilon_t \end{align*}\] where again \(\beta=\alpha \beta^*\) and \(\varepsilon_t \sim \text{NID}(0,\sigma^2)\).

ETS(A,A,N): Specifying the model

ETS(y ~ error("A") + trend("A") + season("N"))

By default, optimal values for \(\beta\) and \(b_0\) are used.

\(\beta\) can be chosen manually in trend().

trend("A", beta = 0.004)
trend("A", beta_range = c(0, 0.1))

Example: Australian population

aus_economy <- global_economy %>% filter(Code == "AUS") %>%
  mutate(Pop = Population / 1e6)
fit <- aus_economy %>%
  model(AAN = ETS(Pop ~ error("A") + trend("A") + season("N")))
report(fit)

Series: Pop 
Model: ETS(A,A,N) 
  Smoothing parameters:
    alpha = 0.9999 
    beta  = 0.3266366 

  Initial states:
     l[0]      b[0]
 10.05414 0.2224818

  sigma^2:  0.0041

      AIC      AICc       BIC 
-76.98569 -75.83184 -66.68347 

Example: Australian population

components(fit) %>% autoplot()

Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_line()`).

Example: Australian population

components(fit) %>%
  left_join(fitted(fit), by = c("Country", ".model", "Year"))

# A dable: 59 x 8 [1Y]
# Key:     Country, .model [1]
# :        Pop = lag(level, 1) + lag(slope, 1) + remainder
   Country   .model  Year   Pop level slope remainder .fitted
   <fct>     <chr>  <dbl> <dbl> <dbl> <dbl>     <dbl>   <dbl>
 1 Australia AAN     1959  NA    10.1 0.222 NA           NA  
 2 Australia AAN     1960  10.3  10.3 0.222 -0.000145    10.3
 3 Australia AAN     1961  10.5  10.5 0.217 -0.0159      10.5
 4 Australia AAN     1962  10.7  10.7 0.231  0.0418      10.7
 5 Australia AAN     1963  11.0  11.0 0.223 -0.0229      11.0
 6 Australia AAN     1964  11.2  11.2 0.221 -0.00641     11.2
 7 Australia AAN     1965  11.4  11.4 0.221 -0.000314    11.4
 8 Australia AAN     1966  11.7  11.7 0.235  0.0418      11.6
 9 Australia AAN     1967  11.8  11.8 0.206 -0.0869      11.9
10 Australia AAN     1968  12.0  12.0 0.208  0.00350     12.0
# ℹ 49 more rows

Example: Australian population

fit %>%
  forecast(h = 10) %>%
  autoplot(aus_economy) +
  labs(y = "Millions", title = "Population: Australia")

Damped trend method

• Damping parameter \(0<\phi<1\).
• If \(\phi=1\), identical to Holt’s linear trend.
• As \(h\rightarrow\infty\), \(\pred{y}{T+h}{T}\rightarrow \ell_T+\phi b_T/(1-\phi)\).
• Short-run forecasts trended, long-run forecasts constant.

Example: Australian population

aus_economy %>%
  model(holt = ETS(Pop ~ error("A") + trend("Ad") + season("N"))) %>%
  forecast(h = 20) %>%
  autoplot(aus_economy)

Example: Australian population

fit <- aus_economy %>%
  filter(Year <= 2010) %>%
  model(
    ses = ETS(Pop ~ error("A") + trend("N") + season("N")),
    holt = ETS(Pop ~ error("A") + trend("A") + season("N")),
    damped = ETS(Pop ~ error("A") + trend("Ad") + season("N"))
  )

tidy(fit)
accuracy(fit)

Example: Australian population

Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing. 
2 observations are missing between 2018 and 2019

term	SES	Linear trend	Damped trend
\(\alpha\)	1.00	1.00	1.00
\(\beta^*\)		0.30	0.40
\(\phi\)			0.98
NA		0.22	0.25
NA	10.28	10.05	10.04
Training RMSE	0.24	0.06	0.07
Test RMSE	1.63	0.15	0.21
Test MASE	6.18	0.55	0.75
Test MAPE	6.09	0.55	0.74
Test MAE	1.45	0.13	0.18

Models with seasonality

Holt-Winters additive method

Holt and Winters extended Holt’s method to capture seasonality.

• \(k=\) integer part of \((h-1)/m\). Ensures estimates from the final year are used for forecasting.
• Parameters:  \(0\le \alpha\le 1\),  \(0\le \beta^*\le 1\),  \(0\le \gamma\le 1-\alpha\)  and \(m=\) period of seasonality (e.g. \(m=4\) for quarterly data).

Holt-Winters additive method

• Seasonal component is usually expressed as \(s_{t} = \gamma^* (y_{t}-\ell_{t})+ (1-\gamma^*)s_{t-m}.\)
• Substitute in for \(\ell_t\): \(s_{t} = \gamma^*(1-\alpha) (y_{t}-\ell_{t-1}-b_{t-1})+ [1-\gamma^*(1-\alpha)]s_{t-m}\)
• We set \(\gamma=\gamma^*(1-\alpha)\).
• The usual parameter restriction is \(0\le\gamma^*\le1\), which translates to \(0\le\gamma\le(1-\alpha)\).

Exponential smoothing: seasonality

ETS(A,A,A)

Holt-Winters additive method with additive errors.

• Forecast errors: \(\varepsilon_{t} = y_t - \hat{y}_{t|t-1}\)
• \(k\) is integer part of \((h-1)/m\).

Holt-Winters multiplicative method

Seasonal variations change in proportion to the level of the series.

• \(k\) is integer part of \((h-1)/m\).
• Additive method: \(s_t\) in absolute terms — within each year \(\sum_i s_i \approx 0\).
• Multiplicative method: \(s_t\) in relative terms — within each year \(\sum_i s_i \approx m\).

ETS(M,A,M)

Holt-Winters multiplicative method with multiplicative errors.

• Forecast errors: \(\varepsilon_{t} = (y_t - \hat{y}_{t|t-1})/\hat{y}_{t|t-1}\)
• \(k\) is integer part of \((h-1)/m\).

Example: Australian holiday tourism

aus_holidays <- tourism %>%
  filter(Purpose == "Holiday") %>%
  summarise(Trips = sum(Trips))
fit <- aus_holidays %>%
  model(
    additive = ETS(Trips ~ error("A") + trend("A") + season("A")),
    multiplicative = ETS(Trips ~ error("M") + trend("A") + season("M"))
  )
fc <- fit %>% forecast()

Example: Australian holiday tourism

fc %>%
  autoplot(aus_holidays, level = NULL) +
  labs(y = "Thousands", title = "Overnight trips")

Estimated components

components(fit)

# A dable: 168 x 7 [1Q]
# Key:     .model [2]
# :        Trips = lag(level, 1) + lag(slope, 1) + lag(season, 4) + remainder
   .model   Quarter  Trips level slope season remainder
   <chr>      <qtr>  <dbl> <dbl> <dbl>  <dbl>     <dbl>
 1 additive 1997 Q1    NA    NA   NA    1512.      NA  
 2 additive 1997 Q2    NA    NA   NA    -290.      NA  
 3 additive 1997 Q3    NA    NA   NA    -684.      NA  
 4 additive 1997 Q4    NA  9899. -37.4  -538.      NA  
 5 additive 1998 Q1 11806. 9964. -24.5  1512.     433. 
 6 additive 1998 Q2  9276. 9851. -35.6  -290.    -374. 
 7 additive 1998 Q3  8642. 9700. -50.2  -684.    -489. 
 8 additive 1998 Q4  9300. 9694. -44.6  -538.     188. 
 9 additive 1999 Q1 11172. 9652. -44.3  1512.      10.7
10 additive 1999 Q2  9608. 9676. -35.6  -290.     290. 
# ℹ 158 more rows

Estimated components

Warning: Removed 4 rows containing missing values or values outside the scale range
(`geom_line()`).
Removed 4 rows containing missing values or values outside the scale range
(`geom_line()`).

Figure 1

Holt-Winters damped method

Often the single most accurate forecasting method for seasonal data:

Holt-Winters with daily data

sth_cross_ped <- pedestrian %>%
  filter(
    Date >= "2016-07-01",
    Sensor == "Southern Cross Station"
  ) %>%
  index_by(Date) %>%
  summarise(Count = sum(Count) / 1000)
sth_cross_ped %>%
  filter(Date <= "2016-07-31") %>%
  model(
    hw = ETS(Count ~ error("M") + trend("Ad") + season("M"))
  ) %>%
  forecast(h = "2 weeks") %>%
  autoplot(sth_cross_ped %>% filter(Date <= "2016-08-14")) +
  labs(
    title = "Daily traffic: Southern Cross",
    y = "Pedestrians ('000)"
  )

Holt-Winters with daily data

Innovations state space models

Exponential smoothing methods

ETS models

Additive error models

Multiplicative error models

Estimating ETS models

• Smoothing parameters \(\alpha\), \(\beta\), \(\gamma\) and \(\phi\), and the initial states \(\ell_0\), \(b_0\), \(s_0,s_{-1},\dots,s_{-m+1}\) are estimated by maximising the “likelihood” = the probability of the data arising from the specified model.
• For models with additive errors equivalent to minimising SSE.
• For models with multiplicative errors, equivalent to minimising SSE.

Innovations state space models

• Estimate parameters \(\bm\theta = (\alpha,\beta,\gamma,\phi)\) and initial states \(\bm{x}_0 = (\ell_0,b_0,s_0,s_{-1},\dots,s_{-m+1})\) by minimizing \(L^*\).

Parameter restrictions

Usual region

• Traditional restrictions in the methods \(0< \alpha,\beta^*,\gamma^*,\phi<1\)(equations interpreted as weighted averages).
• In models we set \(\beta=\alpha\beta^*\) and \(\gamma=(1-\alpha)\gamma^*\).
• Therefore \(0< \alpha <1\),    \(0 < \beta < \alpha\)    and \(0< \gamma < 1-\alpha\).
• \(0.8<\phi<0.98\) — to prevent numerical difficulties.

Admissible region

• To prevent observations in the distant past having a continuing effect on current forecasts.
• Usually (but not always) less restrictive than region.
• For example for ETS(A,N,N): \(0< \alpha <1\) while \(0< \alpha <2\).

Model selection

where \(L\) is the likelihood and \(k\) is the number of parameters initial states estimated in the model.

which is the AIC corrected (for small sample bias).

AIC and cross-validation

Automatic forecasting

From Hyndman et al. (IJF, 2002):

• Apply each model that is appropriate to the data. Optimize parameters and initial values using MLE (or some other criterion).
• Select best method using AICc:
• Produce forecasts using best method.
• Obtain forecast intervals using underlying state space model.

Method performed very well in M3 competition.

Example: National populations

fit <- global_economy %>%
  mutate(Pop = Population / 1e6) %>%
  model(ets = ETS(Pop))

Warning: 1 error encountered for ets
[1] ETS does not support missing values.

fit

# A mable: 263 x 2
# Key:     Country [263]
   Country                      ets
   <fct>                    <model>
 1 Afghanistan         <ETS(A,A,N)>
 2 Albania             <ETS(M,A,N)>
 3 Algeria             <ETS(M,A,N)>
 4 American Samoa      <ETS(M,A,N)>
 5 Andorra             <ETS(M,A,N)>
 6 Angola              <ETS(M,A,N)>
 7 Antigua and Barbuda <ETS(M,A,N)>
 8 Arab World          <ETS(M,A,N)>
 9 Argentina           <ETS(A,A,N)>
10 Armenia             <ETS(M,A,N)>
# ℹ 253 more rows

Example: National populations

fit %>%
  forecast(h = 5)

# A fable: 1,315 x 5 [1Y]
# Key:     Country, .model [263]
   Country  .model  Year             Pop
   <fct>    <chr>  <dbl>          <dist>
 1 Afghani… ets     2018    N(36, 0.012)
 2 Afghani… ets     2019    N(37, 0.059)
 3 Afghani… ets     2020     N(38, 0.16)
 4 Afghani… ets     2021     N(39, 0.35)
 5 Afghani… ets     2022     N(40, 0.64)
 6 Albania  ets     2018 N(2.9, 0.00012)
 7 Albania  ets     2019   N(2.9, 6e-04)
 8 Albania  ets     2020  N(2.9, 0.0017)
 9 Albania  ets     2021  N(2.9, 0.0036)
10 Albania  ets     2022  N(2.9, 0.0066)
# ℹ 1,305 more rows
# ℹ 1 more variable: .mean <dbl>

Example: Australian holiday tourism

holidays <- tourism %>%
  filter(Purpose == "Holiday")
fit <- holidays %>% model(ets = ETS(Trips))
fit

# A mable: 76 x 4
# Key:     Region, State, Purpose [76]
   Region                     State              Purpose          ets
   <chr>                      <chr>              <chr>        <model>
 1 Adelaide                   South Australia    Holiday <ETS(A,N,A)>
 2 Adelaide Hills             South Australia    Holiday <ETS(A,A,N)>
 3 Alice Springs              Northern Territory Holiday <ETS(M,N,A)>
 4 Australia's Coral Coast    Western Australia  Holiday <ETS(M,N,A)>
 5 Australia's Golden Outback Western Australia  Holiday <ETS(M,N,M)>
 6 Australia's North West     Western Australia  Holiday <ETS(A,N,A)>
 7 Australia's South West     Western Australia  Holiday <ETS(M,N,M)>
 8 Ballarat                   Victoria           Holiday <ETS(M,N,A)>
 9 Barkly                     Northern Territory Holiday <ETS(A,N,A)>
10 Barossa                    South Australia    Holiday <ETS(A,N,N)>
# ℹ 66 more rows

Example: Australian holiday tourism

fit %>%
  filter(Region == "Snowy Mountains") %>%
  report()

Series: Trips 
Model: ETS(M,N,A) 
  Smoothing parameters:
    alpha = 0.1571013 
    gamma = 0.0001000991 

  Initial states:
     l[0]      s[0]    s[-1]     s[-2]     s[-3]
 141.6782 -60.95904 130.8567 -42.23776 -27.65986

  sigma^2:  0.0388

     AIC     AICc      BIC 
852.0452 853.6008 868.7194 

Example: Australian holiday tourism

fit %>%
  filter(Region == "Snowy Mountains") %>%
  components(fit)

# A dable: 84 x 9 [1Q]
# Key:     Region, State, Purpose, .model [1]
# :        Trips = (lag(level, 1) + lag(season, 4)) * (1 + remainder)
   Region          State           Purpose .model Quarter Trips level season remainder
   <chr>           <chr>           <chr>   <chr>    <qtr> <dbl> <dbl>  <dbl>     <dbl>
 1 Snowy Mountains New South Wales Holiday ets    1997 Q1  NA     NA   -27.7   NA     
 2 Snowy Mountains New South Wales Holiday ets    1997 Q2  NA     NA   -42.2   NA     
 3 Snowy Mountains New South Wales Holiday ets    1997 Q3  NA     NA   131.    NA     
 4 Snowy Mountains New South Wales Holiday ets    1997 Q4  NA    142.  -61.0   NA     
 5 Snowy Mountains New South Wales Holiday ets    1998 Q1 101.   140.  -27.7   -0.113 
 6 Snowy Mountains New South Wales Holiday ets    1998 Q2 112.   142.  -42.2    0.154 
 7 Snowy Mountains New South Wales Holiday ets    1998 Q3 310.   148.  131.     0.137 
 8 Snowy Mountains New South Wales Holiday ets    1998 Q4  89.8  148.  -61.0    0.0335
 9 Snowy Mountains New South Wales Holiday ets    1999 Q1 112.   147.  -27.7   -0.0687
10 Snowy Mountains New South Wales Holiday ets    1999 Q2 103.   147.  -42.2   -0.0199
# ℹ 74 more rows

Example: Australian holiday tourism

fit %>%
  filter(Region == "Snowy Mountains") %>%
  components(fit) %>%
  autoplot()

Warning: Removed 4 rows containing missing values or values outside the scale range
(`geom_line()`).

Example: Australian holiday tourism

fit %>% forecast()

# A fable: 608 x 7 [1Q]
# Key:     Region, State, Purpose, .model [76]
   Region         State           Purpose .model Quarter
   <chr>          <chr>           <chr>   <chr>    <qtr>
 1 Adelaide       South Australia Holiday ets    2018 Q1
 2 Adelaide       South Australia Holiday ets    2018 Q2
 3 Adelaide       South Australia Holiday ets    2018 Q3
 4 Adelaide       South Australia Holiday ets    2018 Q4
 5 Adelaide       South Australia Holiday ets    2019 Q1
 6 Adelaide       South Australia Holiday ets    2019 Q2
 7 Adelaide       South Australia Holiday ets    2019 Q3
 8 Adelaide       South Australia Holiday ets    2019 Q4
 9 Adelaide Hills South Australia Holiday ets    2018 Q1
10 Adelaide Hills South Australia Holiday ets    2018 Q2
# ℹ 598 more rows
# ℹ 2 more variables: Trips <dist>, .mean <dbl>

Example: Australian holiday tourism

fit %>% forecast() %>%
  filter(Region == "Snowy Mountains") %>%
  autoplot(holidays) +
  labs(y = "Thousands", title = "Overnight trips")

Residuals

Response residuals

\[\hat{e}_t = y_t - \hat{y}_{t|t-1}\]

Innovation residuals

Additive error model: \[\hat\varepsilon_t = y_t - \hat{y}_{t|t-1}\]

Multiplicative error model: \[\hat\varepsilon_t = \frac{y_t - \hat{y}_{t|t-1}}{\hat{y}_{t|t-1}}\]

Example: Australian holiday tourism

aus_holidays <- tourism %>%
  filter(Purpose == "Holiday") %>%
  summarise(Trips = sum(Trips))
fit <- aus_holidays %>%
  model(ets = ETS(Trips)) %>%
  report()

Series: Trips 
Model: ETS(M,N,M) 
  Smoothing parameters:
    alpha = 0.3578226 
    gamma = 0.0009685565 

  Initial states:
     l[0]      s[0]     s[-1]    s[-2]    s[-3]
 9666.501 0.9430367 0.9268433 0.968352 1.161768

  sigma^2:  0.0022

     AIC     AICc      BIC 
1331.372 1332.928 1348.046 

Example: Australian holiday tourism

residuals(fit)
residuals(fit, type = "response")

Example: Australian holiday tourism

fit %>%
  augment()

# A tsibble: 80 x 6 [1Q]
# Key:       .model [1]
   .model Quarter  Trips .fitted .resid   .innov
   <chr>    <qtr>  <dbl>   <dbl>  <dbl>    <dbl>
 1 ets    1998 Q1 11806.  11230.  576.   0.0513 
 2 ets    1998 Q2  9276.   9532. -257.  -0.0269 
 3 ets    1998 Q3  8642.   9036. -393.  -0.0435 
 4 ets    1998 Q4  9300.   9050.  249.   0.0275 
 5 ets    1999 Q1 11172.  11260.  -88.0 -0.00781
 6 ets    1999 Q2  9608.   9358.  249.   0.0266 
 7 ets    1999 Q3  8914.   9042. -129.  -0.0142 
 8 ets    1999 Q4  9026.   9154. -129.  -0.0140 
 9 ets    2000 Q1 11071.  11221. -150.  -0.0134 
10 ets    2000 Q2  9196.   9308. -111.  -0.0120 
# ℹ 70 more rows

Some unstable models

• Some of the combinations of (Error, Trend, Seasonal) can lead to numerical difficulties; see equations with division by a state.
• These are: ETS(A,N,M), ETS(A,A,M), ETS(A,A,M).
• Models with multiplicative errors are useful for strictly positive data, but are not numerically stable with data containing zeros or negative values. In that case only the six fully additive models will be applied.

Exponential smoothing models

Forecasting with exponential smoothing

Forecasting with ETS models

iterate the equations for \(t=T+1,T+2,\dots,T+h\) and set all \(\varepsilon_t=0\) for \(t>T\).

• Not the same as \(\text{E}(y_{t+h} | \bm{x}_t)\) unless seasonality is additive.
• fable uses \(\text{E}(y_{t+h} | \bm{x}_t)\).
• Point forecasts for ETS(A,*,*) are identical to ETS(M,*,*) if the parameters are the same.

Forecasting with ETS models

can only be generated using the models.

• The prediction intervals will differ between models with additive and multiplicative errors.
• Exact formulae for some models.
• More general to simulate future sample paths, conditional on the last estimate of the states, and to obtain prediction intervals from the percentiles of these simulated future paths.

Prediction intervals

Example: Corticosteroid drug sales

h02 <- PBS %>%
  filter(ATC2 == "H02") %>%
  summarise(Cost = sum(Cost))
h02 %>% autoplot(Cost)

Example: Corticosteroid drug sales

h02 %>%
  model(ETS(Cost)) %>%
  report()

Series: Cost 
Model: ETS(M,Ad,M) 
  Smoothing parameters:
    alpha = 0.3071016 
    beta  = 0.0001006793 
    gamma = 0.0001007181 
    phi   = 0.977528 

  Initial states:
     l[0]    b[0]      s[0]     s[-1]     s[-2]     s[-3]     s[-4]    s[-5]    s[-6]
 417268.7 8205.82 0.8716807 0.8259747 0.7562808 0.7733338 0.6872373 1.283821 1.324616
    s[-7]    s[-8]    s[-9]   s[-10]    s[-11]
 1.180067 1.163601 1.104801 1.047963 0.9806235

  sigma^2:  0.0046

     AIC     AICc      BIC 
5515.212 5518.909 5574.938 

Example: Corticosteroid drug sales

h02 %>%
  model(ETS(Cost ~ error("A") + trend("A") + season("A"))) %>%
  report()

Series: Cost 
Model: ETS(A,A,A) 
  Smoothing parameters:
    alpha = 0.1702163 
    beta  = 0.006310854 
    gamma = 0.4545987 

  Initial states:
     l[0]     b[0]      s[0]     s[-1]     s[-2]     s[-3]     s[-4]    s[-5]    s[-6]
 409705.9 9097.111 -99075.37 -136602.3 -191496.1 -174530.8 -241436.7 210643.8 244644.2
    s[-7]    s[-8]    s[-9]  s[-10]    s[-11]
 145368.2 130569.6 84457.69 39131.7 -11673.71

  sigma^2:  3498869384

     AIC     AICc      BIC 
5585.278 5588.568 5641.686 

Example: Corticosteroid drug sales

h02 %>%
  model(ETS(Cost)) %>%
  forecast() %>%
  autoplot(h02)

Example: Corticosteroid drug sales

h02 %>%
  model(
    auto = ETS(Cost),
    AAA = ETS(Cost ~ error("A") + trend("A") + season("A"))
  ) %>%
  accuracy()

Model	MAE	RMSE	MAPE	MASE	RMSSE
auto	38649.04	51102.24	4.988983	0.6375806	0.6891173
AAA	43378.40	56784.23	6.047574	0.7155993	0.7657394